c   MEPU.for generates uncertainty estimates

c    at 95, 50, 5th fractiles and the corresponding size

c    distribution.  Note, distribution is average of size at

c    fractile +/- 10 observations on either size.

c    These uncertainty estimates are for recoverable oil

c    or NA gas for ANWR individual plays. Program could be easily

c    modified to get in-place estimates of uncertainty.

c

c    Written by Schuenemeyer 3/22/98

c

c   Input files:

c    Unit  8 - play file

c    Unit  9 - prospect file

c   Output files:

c    Unit 10 - uncertainty estimates

c     Note: output will be pasted into ANWR Summary & Distns

c      Excel spreadsheets

c

c   Subroutines required:

c     SizeClass, Buble, Bublei(included in this program)

c

      character*25 fnpl,fnpr,fndn

	character title*80

      dimension x(10000),id(10000),tmp(7),nar(30),naid(30),so(2,13)

	dimension nt(4),sm(2),perc(3)

	integer fsc

	data perc/0.05,.5,.95/

	call getdat(iyrx,imonx,idayx)

	marry=13

	nsper=10

      write(*,*)' Enter play file name'

	read(*,'(a25)')fnpl

      open(8,file=fnpl)

	write(*,*)' Enter prospect file name'

	read(*,'(a25)')fnpr

      open(9,file=fnpr)

	write(*,*)' Enter output uncertainty file name'

	read(*,'(a25)')fndn

      open(10,file=fndn)

    7	write(*,8)

    8 format(/' Enter total number of PLAYS & Oil/Gas(1or2)')

	read(*,*) npl,nog

	if(nog.lt.1.or.nog.gt.2) goto 7

	if(nog.eq.1) then 

	 scdiv=1.

	 idv=1

	  else

	 scdiv=6.

	 idv=3

	end if

	knt=0

c  read oil or gas values from play file

	read(8,*)

   10	read(8,*,end=20)it,(tmp(i),i=1,7)

	knt=knt+1

	id(knt)=it

	x(knt)=tmp(idv)

	goto 10

   20	close(8)

      iexc=npl-knt

	if(iexc.lt.0) pause 1

	nbs=0

c  big percentile (fractile) loop

      do ki=1,3

      ndep=0

      do i=1,2

	 sm(i)=0

	 do j=1,marry

	  so(i,j)=0.

	 end do

	end do

	maxf=0

      if(nbs.eq.0) then

c   sort oil/gas play totals

	 call buble(x,id,knt)

	nbs=1

	end if

      npid=perc(ki)*npl+.5

c   percentile is zero

	if(npid.le.iexc) goto 56

c   relate npid to actual data and interval bounds

	ndat=npid-iexc

	une=x(ndat)

	 write(*,*) ' Uncert est',une

c   get size distribution

      nll=ndat-nsper

	if(nll.lt.1)nll=1

	nul=ndat+nsper

	if(nul.gt.knt)nll=knt

	ndif=nul-nll+1

	xndif=ndif

c   store id's

      do i=nll,nul

	 ia=i-nll+1

	 nar(ia)=id(i)

	 naid(ia)=ia

c	write(*,*)ki,i,ia,nar(ia)

	end do

      call bublei(nar,naid,ndif)

c   nar contains the sorted number for the prospect file.

c   Get play data and summarize.

c

c   read headers from prospect file

	rewind 9

      read(9,'(a70)')title

	read(9,*)nvar,nobs

	do i=1,nvar

	 read(9,*)

	end do

	nsr=0

      do i = 1,ndif

	 ia=nar(i)

	if(nsr.eq.1) goto 30

   28	read(9,*,end=59) (nt(k),k=1,4),oilg

   30 nsr=0

	if(nt(1).lt.ia) goto 28

	if(nt(1).gt.ia) then

	 nsr=1

	 goto 33

	end if

	if(nt(3).eq.1) nt2=nt(2)

	if(nt(4).eq. nog) then 

	pet=oilg/scdiv

	call SizeClass(pet,maxf,marry,fsc)

c  accumulate

	ndep=ndep+1

      so(1,fsc)=so(1,fsc)+1

	so(2,fsc)=so(2,fsc)+oilg

	end if

	nt2=nt2-1

	if(nt2.gt.0)goto 28

   33	end do

c  summarize and write out results

	do j=1,marry

	do i=1,2

	 so(i,j)=so(i,j)/xndif

	 sm(i)=sm(i)+so(i,j)

	end do

	end do

   56	fra=(1.-perc(ki))*100.

      write(10,60)fra

   60 format(f15.0)

   61	do i=1,2

	write(10,65)sm(i),(so(i,j),j=1,marry)

   65 format(3x,14f12.3)

      end do

      write(*,*) ' Max size class',maxf

	end do

      close(10)

	stop

   59 pause 3

      END



	Subroutine SizeClass(size, maxf, marray, fsc)

c Computes size class in log based 2

      integer fsc

       fsc = Int(Log(size) / Log(2.)) - 1

       If (fsc .lt. 1) fsc = 1

       If (fsc .gt. maxf)  maxf = fsc

       If (fsc .gt. marray)  fsc = marray

	return

      End 



	SUBROUTINE BUBLE(X,ID,N)

      DIMENSION X(1),ID(1)

      KS=N

   15 KW=0

      DO 30 I=2,KS

      IF(X(I).GE.X(I-1)) GOTO 30

      TMP=X(I)

      X(I)=X(I-1)

      X(I-1)=TMP

      NTI=ID(I)

      ID(I)=ID(I-1)

      ID(I-1)=NTI

      KW=1

   30 CONTINUE

      IF(KW.EQ.0) RETURN

      KS=KS-1

      IF(KS.EQ.1) RETURN

      GOTO 15

      END



	SUBROUTINE BUBLEI(X,ID,N)

      integer x,tmp

      DIMENSION X(1),ID(1)

      KS=N

   15 KW=0

      DO 30 I=2,KS

      IF(X(I).GE.X(I-1)) GOTO 30

      TMP=X(I)

      X(I)=X(I-1)

      X(I-1)=TMP

      NTI=ID(I)

      ID(I)=ID(I-1)

      ID(I-1)=NTI

      KW=1

   30 CONTINUE

      IF(KW.EQ.0) RETURN

      KS=KS-1

      IF(KS.EQ.1) RETURN

      GOTO 15

      END

